Biphonons in the Klein-Gordon lattice 



Laurent Proville 

Service de Recherches de Metallurgie Physique, CEA/DEN/DMN Saclay 91191-Gif-sur-Yvette Cedex, France 

(Dated: February 2, 2008) 

A numerical approach is proposed to compute the phonon bound states in a quantum nonlinear 
Klein-Gordon lattice. In agreement with other studies»*-on a different quantum lattice, nonlinearity 
is found to lead to a phonon pairing and consequently some biphonon excitations. The energy branch 
and the correlation properties of the Klein-Gordon biphonon are studied in detail. 

PACS numbers: 63.20.Ry, 03.65.Ge, ll.10.Lm, 63.20.Dj 

I. INTRODUCTION 

In lattices made of identical particles, the energy is formulated by the Hamiltonian operator: 

l j=<l> 

where xi and pi are displacement and momentum of the particle at site I, in a d-dimensional lattice. From the left 
to the right hand side of the equation Eq^ the energy contributions are identified as the kinetic energy, the local 
potential and the interaction between particles. Our purpose is to study the case of quantum particles that are weakly 
interacting, i.e., the onsite energy V dominates the interaction W. Physically, this can account for the coupling of 
internal modes in molecular crystals. For small amplitudes of xi, the well-known harmonic approximation reduces 
H to a sum of quadratic terms, i.e., the linear Klein-Gordon (KG) Hamiltonian. So, the Schrodinger equation can 
be solved analytically. The elementary excitation is a plane wave called an optical phonon and whose energy is fixed 
by the wave momentum q, in the lattice Brillouin zone. A consequence of the ideal harmonicity is that higher order 
excitations are simply the linear superpositions of these optical phonons. 

For larger displacements, some non-quadratic contributions are involved in the expansion of H. Then the nonlinear 
KG Hamiltonian can no longer be diagonalized analytically. In a nonlinear KG lattice, F. Bogani^ derived some one- 
and two- phonon renormalized Green functions and showed that the nonlinear terms involve a pairing of the optical 
phonon modes. It confirmed the existence of biphonon excitations, studied earlier in a different lattice model by V.M. 
AgranovicfA A convincing agreement was found between theorjji and experiments^ 6 in molecular crystals where the 
internal molecule bonds yield a strong nonlinearity. The direct diagonalization of a KG Hamiltonian is, in principle, 
more precise than the computation of Green functions since it requires fewer approximations. In Ref. [7], by treating 
numerically the KG model, W. Z. Wang and al. confirmed the existence of phonon bound states. Furthermore, some 
of these states have been shown to feature a particle-like energy band, for certain model pa rameters. The authors 
identified these specific excitations as being some quantum breathers (see Refs. I2u8ll9ll0llll for more details about 
quantum breathers) because of their counterparts in classical mechanicsiS* 1 ^. Nonetheless, the approach proposed in 
Ref.Q requires a huge computing cost so the size of lattices was limited to a one-dimensional (ID) chain of 8 unit 
cells. Moreover, the numerical simulations were restricted to the parameter region where the non-harmonic part of 
the lattice energy is modelled by a quartic onsite potential, i.e., the well-known 4 model. In the present paper, we 
propose a numerical treatment of the nonlinear KG lattice which takes advantage of the weak coupling [W in Eq^l . 
That permits to analyze lattices large enough to approach the infinite system features and to study different types of 
nonlinearity, as well as the two-dimensional (2D) KG lattice. 

We confirm that when nonlinearity is significant, a pairing of optical phonon states occurs and the so-called 
biphononi branch contributes to the energy-spectrum. That branch splits from the two-phonon band by opening 
a gap. The width of that gap indicates the magnitude of nonlinearity since the biphonon gap vanishes completely 
for a pure harmonic lattice. In between the two types of lattice, i.e., harmonic and strongly nonlinear, the binding 
energy of the biphonon drops to zero at the center of the lattice Brillouin zone (BZ) while at the edge, the biphonon 
excitations are still bound. Then, in the energy-spectrum, the biphonon gap vanishes at the center, whereas a pseu- 
dogap is found to open at the edge of BZ. We predict that the pseudogap is a systematic feature of lattices in which 
nonlinearity is moderate, whatever is the lattice dimension or the dominant nonlinearity, i.e., </> 3 or 4 . In addition, 
we enhance how quantum properties of the biphonon depends on the nonlinearity. When the biphonon gap opens, 
the Klein-Gordon biphonon excitations show a finite correlation length, for all momentum q, under the condition that 
the non-quadratic energy term is a 4 potential. That agrees with findings of Ref 7. Considering the cubic term in 
the potential energy V, it involves a long range correlation of the biphonon states. The space correlation properties 
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of triphonon are also studied and our results are used to establish some expectations on the existence of breather-like 
excitations in the quantum KG lattice. 

The present paper is organized as follows. In Sec[H]the model for the nonlinear discrete lattice is introduced and 
the computing method is detailed. Then it is tested for the quadratic lattices as well as for the 4> 4 model. In Sec lIIII 
our results are presented concerning the biphonon spectrum while in Sec II VI the space correlation properties of the 
phonon bound states are studied. Finally, these results are discussed in SecfVl 



II. MODEL AND COMPUTING METHOD 



In Eq^ at node I of a translational invariant d-dimensional lattice, the quantum particle of mass m evolves in a 
local potential V, being coupled to its nearest neighbours, j by the interaction W. For moderate amplitudes of mass 
displacements around equilibrium, V and W can be expanded as Taylor series. The expansion of V is truncated 
to the fourth order V(xi) — a,2X 2 + a^xf + a^xf while for W, only the quadratic term is retained, W(xi — Xj) — 
—c(xi — Xj) 2 . Higher order terms can be treated with no difficulty in what follows. Actually, they are found not to 
change qualitatively the results, at least for reasonable val ues of energy coeffici ents, co nsistent with optical modes. 

Introducing the dimensionless operators Pi = pi/ %/ mhfl, Xi = xi yj mfl/h where the frequency Q = 
\j2{a,2 — 2.c.d)/m is defined for either a chain d = 1 or a square lattice d = 2, the Hamiltonian reads 

H = hnY J ^ + ^+A,Xf + A i Xt + ^Xi X i (2) 
i ]=<i> 

where one finds the dimensionless coefficients: 



4c 

^3 = fl3\/ — A± = a 4 — — - and C = 



i 3 n 5 ' * m 2 n 3 mil 2 ' 

For the harmonic lattice, i.e., A3 = and A4 = 0, the Fourier transform of both the displacements, Xi = 
-^J2 q e ~ l9Xl Xq and the momenta, Pi = -^J2 q e ~ lqxl P q simplifies H into a sum of independent Hamiltonian: 

h q = hSl(^ + ^lx 2 ) (3) 

with uj q = \/T— 2C^2 k cos(qk) and qk is the dimensionless coordinate of the wave vector in the k th direction of the 
lattice. The periodic boundary conditions impose qk — 2irlk/ Lk where is the number of sites in the fc th direction 
and lk is an integer lk € [1,-kfc]- The lattice size is denoted S = Hf =1 Lk- Using the standard harmonic oscillator 
theory, one finds the eigenvalues: 




A {nJ = fto^K+9)W 1 + 2C E cos (*)- ( 4 ) 



The n q are quantum numbers that range from to infinity and fix the energy contribution of the mode q. For 
the nonlinear case, the aforementioned procedure is no longer simple because non-quadratic terms yield a coupling 
between the h q operators. Indeed, writing the nonlinear energy term Q^j Xf) as a function of the displacements 
Fourier transform X q , gives (-7= J2 q q > X q X q > X^_ q _ ql )). The computation of the corresponding bracket thus scales as 

S 2 as the sum runs over 2 wave vector, while for the quartic term it would scale as S 3 . Such a task has been achieved 
in Ref|3 for the quartic term. In the algorithm which follows, the computation of the brackets in Eq|Sl scales as S 
which requires much less computation time for a given S. Starting from the exact diagonalization of the Hamiltonian 
where no interaction couples displacements, the low energy states are used to construct a set of Bloch waves upon 
which is expanded the entire Hamiltonian, including the coupling W. The Schrodinger equation is then approximately 
solved with an error which shrinks to zero by increasing the basis cutoff. 

pi . 

The starting point is thus the eigenvalue problem for a sing le oscillator hi = + ^+ A 3 Xf + AaX?). The 

Bose-Einstein operators af = (A; — iP;)/v2 and a; = (A; + iP;)/V2 are introduced in the writing of hi: 



1 A-i 

hn[a+ai + - + -4=(af 3 + a f + 3a+af + 3a+ 2 a ; + 3a+ + 3a z ) 

2 y8 

A-4 / +4 1 4 , 
— ( 0,, -I- a.T 4- 
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+ '-iriat 4 + a t + W 3a i + W a f + 6 ( a /" 2a ? + a F* + a f + 2a t a l) + 3 )1- ( 5 ) 
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Expanding the operator hi on the Einstein states, i.e., \n, I >= -h=<2; \$i > for all n £ {0...N — 1} gives a matrix M. 
of rank N. In each row of A4, one finds the nonzero coefficients: 

3 3 1 

M n , n = -A 4 n 2 + n{3A 4 + 1 ) + -A 4 + - 

M n ,n+4 = -At V(" + 4)(n + 3)(n + 2)(n + 1) 

■M n , n+3 = -^=v4 3V /(n + 3)(n + 2)(n + l) 
3 



M„, n +i = -j=A^(n + If 



M n , n+2 = -Ai(An + 6)y/(n + 2)(n + l). (6) 

When the cutoff N tends to infinity, the Einstein states form a basis in the space of onsite states. The Schrodinger 
problem for the Hamiltonian hi is thus equivalent to the diagonalization of M.. That diagonalization is compared to 
the semi-classical quantization^ in Fig. ^a), for the case of a He atom embedded into a double- well potential: V(x) = 
\&Ed/b i x 2 {x — b) 2 . The parameters, Ed and 6, are the energy barrier and distance between minima, respectively (see 
Ref. Ila) . The very good agreement proves the efficiency of the diagonalization method even for a non-monotonic 
onsite potential. Arranging the onsite eigenstates in increasing order of their eigenvalues, the a th eigenstate is denoted 
4> a ^i and its eigenvalue is 7(a). As shown in Fig^b), each eigenvalue 7(a) is found to converge to a steady value as N 
increases. In Fig^(b), the graphic does not allow to distinguish N = 00 from N > 100. With today's computers, the 
cutoff has been easily increased to N = 2000 which has been taken for the limit N = 00 in Fig^b). Increasing the 
cutoff to values larger than N — 100 does not change significantly our final results on the low energy excitations of the 
KG lattice fSec lIIIfl . In what follows, N is cautiously fixed to N = 500 and then the time requirement to diagonalize 
A4 is about few minutes on a PC computer. Note that increasing N implies no overload of the calculations in the 
second part of the algorithm. 

We now treat lattices with non-zero inter-site coupling. The onsite state products Hi4> ai ,i form a complete orthogonal 
base for the lattice states. In order to reduce the computer memory requirement, one takes advantage of the translation 
invariance by introducing the Bloch wave formulation for the state products. Among those states some equivalence 
classes can be constructed in which each state results from a translation applied to another state of the same class. 
Retaining only one element for each translation class, the state which represents the class is identified by the series of 
its aj's, that is denoted [Iljaj]. The construction of the equivalence classes is performed numerically. For each class, 
a Bloch wave can be written as follows: 

B [Uiat] (q) = * Yl e~ iM 'ni4w +j (7) 

where Arn iai i ensures the normalization. Some attention must be paid to the possible translation symmetry of the state 
products that may be higher than the lattice symmetry. Indeed, for a given product there may be a lattice vector t that 
verifies Hi<f> ai ,i — ^i^ai.i-t with coordinates tk such as t% < Lp.. It implies that Am iai ] = S.IVf. =1 {Lk/tk — fc(Lk/tk)) 
where fc(Lk/tk) is the fractional portion of the ratio L^/tk- Then the Bloch wave can only take the momentum q 
such as qk = 2irpk/Lk = 2irp' k /tk where pk and p' k are some different integers. The set of states {Bm ia A (q)} q ,N cut , 
including the uniform state Ilj^o.i at q = Q, form a truncated basis where N cut fixes the upper boundary on the onsite 
excitations: J2i a i — N cut . When C is negligible, these states are the eigenstates of H. For moderate values of C, 
they should be good approximates. Since the Bloch waves with different q, are not hybridized by H, the Hamiltonian 
can be expanded separately for each q. It is performed analytically and gives a matrix B{q) the coefficients of which 
are written as follows: 

< B [nzai] (q)\H\B {I] _ z0z] {q) >= =[H t (? ai A V 7(0^) - — Vexp(-ig x j) 

^A [n . ai] A [n .p i] ; I - 



D(ai, 0i + j)D(o4 + k,0i+k+j)'Hi?i,i+ktiai,Pi+ j ] (8) 



k=<l> 

where D(ai,@i) denotes the bracket < (fra^X^^ > that is given by: 

1 N 

Dfafo) = -t=Yj < M 1 ' i > + < 1 + h*\4>p,i > <l~h #/3,« >)■ (9) 

v 2 1=0 
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The eigenvalues of B(q) are computed numerically with an exact HouseHolder method 16 . In Fig. [3 for a ID chain 
S = 4, our calculation is compared to Ref. (see Ref . Il7l for conversion of model parameters). A very good agreement 
is noted for the low energy states since the eigen-spectra are superposed in Fig.0 In Ref.0, the Schrodinger equation 
was solved by diagonalization of the matrix obtained from expanding H in the Einstein phonon basis, i.e., the 
eigenstates of the pure harmonic lattice (see Eq|3J). According to the authors : "it restricts the numerical simulations 
to a parameter region where nonlinearity is not too large". In contrast, thanks to the first step of our algorithm which 
solves the single site nonlinear eigenvalue problem, we can treat all types of nonlinearity (weak or strong and with 
<f> 3 or tfi 4 terms), provided the inter-site coupling is not too large. For instance, the approach of Ref|7] requires some 
computations even for C = if ^3 7^ or A4 ^ 0, which is straightforwardly solved in the first step of the present 
algorithm. Within the second step, in order to estimate the accuracy of our calculations, different lattice sizes have 
been tested for a reasonable value of C (C — 0.05/d, see Sec. IIII|I and different A 3 and A4. For small sizes, N cut can 
be stepped up sufficiently to make eigenenergies converge to steady values. The convergence is as fast as C is small, 
i.e., when C = the computation is exact (to machine precision) and near instantaneous whereas when C is raised, 
the accuracy becomes worse because of inter-site coupling terms: (afa'j) and (a^a?) that involve hybridization with 
high energy Bloch waves (Eq[7|), above the cutoff. Once N cut has been determined to achieve the required precision, 
then the lattice size is increased up to the capacity of our computer memory. For instance, in a ID lattice with 
S = 17, N cut has been varied from N cu t — 3 (68 Bloch waves, EqEJ) to N cut — 6 (5940 Bloch waves). For N cut — 4, 
the error on the low energy eigenvalues, say the 2 phonon states, is inferior to 1% in comparison with N cut = 6. The 
size has thus been increased to S = 33 (3052 Bloch waves) with no noticeable discrepancy of the eigen-spectrum. For 
such a lattice, the time required for the computation of the matrix B{q) scales in minutes whereas the diagonalization 
requires few hours with a PC. This can be reduced to some minutes with a vectorial computer and suitable numerical 
libraries. The matrices we have to treat are much smaller than those in Ref 7, which accounts for the tractability of our 
method. For 5 = 4, 65536 states were required in Ref0 whereas only 19 Bloch waves were required in our calculations 
for the same lattice, with same parameters (see Fig|2J). This increase in efficiency has been possible because we took 
advantage of the weak inter-site coupling. For the 2D lattices, the size has been limited to S = 13 x 13 which involves 
4931 Bloch waves with N cut = 3. Some improvements are under investigation. For example, the number of required 
states can be reduced again by imposing a* < ni ow for Bloch waves that verify J^i a i > Ni ow in Eq0 (A/ oto < N cut ). 
The integers ni ow and Ni ow are then adjusted so as not to change the precision over the low energy eigenstates. 



III. GAP AND PSEUDOGAP IN THE OPTICAL PHONON SPECTRUM 



For different values of nonlinear coefficients A3 and A4, we first examine the vibration spectrum of a one-dimensional 
(ID) lattice. The 2D lattice is treated at the end of the present section. When the non-quadratic part of the lattice 
energy is negligible, the eigen-spectrum of H is composed of the fundamental optical branch due to the harmonic 
phonon states (in Eq. a single q verifies n q = 1) and the branches due to the linear superpositions of these phonons 
(in Eq. 0] several g's verify n q = 1). The latest branches are stacked together into distinct bundles, each of them 
filling in a compact range of energy. In Fig. [3] and following ones, each eigenvalue of the finite size Hamiltonian is 
plotted as a single circle symbol. The distinct eigen-energies participate in different branches. The phonon branch is 
marked with the tag {1} while the branches that are due to the linear superposition of 2 phonon states are labelled 
by the tag {11}. For a macroscopic system, the bundle {11} covers a dense range of energy and forms a continuous 
band. The width of an optical phonon branch being physically a few percent of the elementary excitation energy, the 
dimensionless coupling C is fixed to C = 0.05 which gives, indeed a phonon branch width Ai 10% of the phonon 
energy (see left inserts in Figs. I3I4|I . When nonlinearity is significant, the phonon branch shows no qualitative change 
(left inserts in Figs-OJa) and^Ja)) in comparison with the fundamental optical branch in harmonic lattice. On the 
other hand, an isolated spectral branch is found in addition to the phonon branch and its combination tones (see 
FigsGJa) and^Ja) and right-hand inserts). In Fig[5l varying artificially the coupling C from the trivial case C = 0, 
demonstrates that the additional branch coincides with the energy of the Bloch wave Br a o,...,o] w ith a single onsite 
excitation a = 2. In Fig|3]and the next ones, the additional branch is marked with a single tag {2}. By analogy with 
biphonon theory!, this branch is identified as the biphonon energy. Similar results are found for the triphonon states 
whose branch is labelled by the tag {3} in Fig[S] The reason for these isolated branches is that onsite Hamiltonian 
eigenvalues r y a >i do not match the linear fit given by (71 — "/o)a + 70. This is the consequence of hi anharmonicity. 
The differences 7 Q >i — [(71 — 7o)a + 70] involve some gaps in the C = spectrum which is composed of the Bloch 
wave energies. A moderate inter-site coupling hybridizes these states Eq[7|but the largest gaps remain (Fig. [SJ. The 
raising of degeneracy of B\n iai ] where only 2 a^s equal 1 and the rest are zero (number of these states is S(S — l)/2, 
their energy is 27(1) + (S — 2)7(0) at C — 0) yields a bundle of branches which correspond to the linear superposition 
of 2 single phonon states. In Fig0 for higher energy, other bundles are labelled out by the tags {111} and {21}. 
At zero coupling, these branches coincide with the energies of states -B[i.i,i.o,...,o] an d -Bp, 1,0,.. .,()]■ For a macroscopic 



5 



lattice, the bundles {111} and {21} form some dense bands, as well as {11}. They are the unbound associations of 3 
phonon states and of a biphonon with a single phonon, respectively. In FigsEJa) and 01 a), for different parameters, 
the biphonon branch splits from the 2 phonons band. Measuring the energy of a biphonon state with reference to the 
unbound 2 phonons for same momentum q, a binding energy of biphonon is defined. A positive binding energy occurs 
when the onsite potential V is harder than a harmonic function (Fig|3Ja)) whereas a softening yields a negative binding 
energy (FigQJa)). The biphonon energy gap is determined as the minimum of the absolute value of binding energy 
with respect to q. The biphonon gap reveals the strength of nonlinearity since when the biphonon gap overpasses the 
phonon branch width (as it does in Figs. |3J a) and^Ja)), it clearly indicates a significant contribution of non-quadratic 
terms. While in FigslSfa) and^Ja), a biphonon gap opens, it is found that when nonlinearity is weak the biphonon 
binding energy vanishes at center of the lattice Brillouin zone (BZ) (FigsOJb) and^Jb)). However, at the edge of BZ, 
the binding energy is comparable to the width of the phonon branch Ai (inserts in Figs|3fb) and^Jb)). Consequently, 
a pseudogap is yielded when the non-quadratic energy has same magnitude as inter-site coupling. In this regime, the 
biphonon excitations exist only at the edge of BZ while they are dissociated into unbound phonons at center. With 
similar results, other calculations have been performed for different parameters. They showed that the biphonon 
pseudogap is a generic feature of lattices where nonlinearity is moderate. Similar pseudogaps have been noted in 
different quantum Iatticesi2i22*21. The pseudogap opens at the edge of BZ, even though the coupling sign is changed. 
So it is the g-range where nonlinear behavior is likely to be experimentally measured in materials where nonlinearity 
is weak. 

In Fig|31 the low energy eigenvalues of H are plotted versus parameter C. The variations of the energy branches 
of the nonlinear excitations are labelled by tags defined previously. It can be noted that the widths of bound states 
branches increase with C much slower than unbound phonon bands. The branch of the a phonon bound states 
(biphonon for a — 2 and triphonon for a = 3) are found to merge with unbound phonon bands for a certain threshold 
C a . At C < C a , the a th branch and the unbound phonon bands are separated by a gap whereas around C ~ C a , 
only a pseudogap separates them partially The C a threshold depends on both coefficients A 3 and A 4 and it is 
different for each a phonon bound state branch because of anharmonicity. A unique set of nonlinear parameters 
A% and A4 corresponds to the energy distribution of the biphonon and triphonon branches. So in principle, if a 
spectroscopy is able to measure the biphonon and the triphonon resonances, it is sufficient for inverting our numerical 
treatment and thus determining the nonlinear parameters. Moreover, the results in Fig|5] demonstrate that even 
though the non-quadratic terms in V are not large enough to open a biphonon gap, i.e., C > C2 then a gap or at 
least a pseudogap opens for the a phonon bound states with a > 2. Finally, the theoretical results in Fig|3Jb) are 
qualitatively similar to the experimental findings in Refll8l in which the Raman analysis of a molecular H2 crystal 
shows a pressure-induced bound-unbound transition of the so called bivibron around 25 GPa. There is indeed a 
likeness between Fig. 10 in Ref . Il8l and the 2 phonon energy region in Fig.[SJb). In our model, the pressure variation of 
experiments^ can be simulated by a change of the coupling parameter C due to the fact that neighbouring molecules 
are moved closer together because of the external pressure. Actually, the increase of C induces a bound-unbound 
transition of the biphonon at C — Ci- 

In Fig- El the diagonalization of a 2-dimensional (2D) lattice Hamiltonian is performed for A3 = and for different 
values of A4. The coupling amplitude is such as the phonon band width Ai is a few percent of the elementary excitation 
branch. By estimating that Ai w 2.<£|C|, the value of the dimensionless coupling is fixed around C = 0.025. In the 
first overtone region, a gap opens when A4 is large (top of FigEJ) whereas that gap closes at the center of the BZ 
when A 4 — 0.025 (bottom of Fig. EJ). In the latter case, a pseudogap is found to open around q = [11] and the width 
of that pseudogap has same order of magnitude as the phonon band width. These both quantities can be compared 
in the inserts of Fig|S] where the spectrum profile along [11] is plotted. The pseudogap width is same as in Fig|3]for 
the ID chain. Consequently, a pseudogap is expected for all lattice dimensions when both the inter-site coupling and 
the non-quadratic energy have a comparable magnitude. 

IV. CORRELATION PROPERTIES OF THE PHONON BOUND STATES 

In the present section, the study is focused on the space correlation in a ID lattice. The space correlation function 
for displacements is defined as follows: 

/(*, n) = J2 < *|^f+n|# > - < $l*z|$ >< *|Jfi+n|* > (10) 
1 

where $ is an eigenstate and n is the dimensionless distance (n > 0). For an harmonic lattice, at weak inter-site 
coupling, it is easily verified that the single phonon correlation function is well approximated by [cos{q n)). Then, 
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the correlation functions exhibit a space modulation with a constant amplitude. At q — 0, the function /($,ra) is a 
non-zero constant for all ra. 

Before discussing our results, it is convenient to compute analytically the function /(<&, ra) for the Bloch waves 
B[ Uiai ] (Eq. 0) with only one onsite excitation of order a. These states have the form: 

B[ a .o,...fi](q) = -j^^e'^^k+j.Ili^kfai+j. (11) 

3 

Let us introduce the notation, ipaiq) for these states. At the uncoupled limit (C = 0), the onsite excitation Bloch 
waves correspond to some eigenstates of H (see Fig|S]and text in Sec lIII|) . The space correlation function of ipaiq) is 
given by: 

/CM*), n) = 2D(a, 0)» cos ( 9 ra) - ^(a,a)-D { 0^ (12) 

where D(a,0) and D(a,a) are defined in EqEl The constant term in the right hand side of Eq. 1121 shrinks to zero 
when the lattice size S increases. In the numerical computations that are presented below, for finite lattices, that 
constant is removed in order to approach the behavior of the infinite lattices. Another point which is noteworthy 
is that, in general, the displacement operators are correlated for a given state i/) a (q). Indeed the bracket D(a,0) is 
usually not zero, except for even values of a when A3 = 0, and the function f(tpa(l), n) shows a space modulation of 
amplitude R a = 2D(a,0) 2 . The coefficients R a are given in table Tab[I]for parameters of Figs|3and2J 

The function /(<&, ra) is plotted in Figs 17191 for eigenstates $ that are either the phonon, biphonon or the triphonon 
states. For each of these states, /($,ra) is computed for several values of q. The calculations at C > 0, demonstrate 
that the long range behavior of the correlation function for the a phonon bound states (1 < a < 4) as well as for 
the phonon states (a — 1) is similar to the variations of the function f(ip a (q),n), i.e., the modulations have same 
amplitudes for large distances. For a weak coupling, this property holds provided that the a th discrete branch is 
separated from the rest of energy spectrum by some gaps. In Figs[7{a)-02a), the correlation function of phonon states 
shows a modulation which extends over the whole lattice as in the pure harmonic lattice. For large ra, the variations 
of /($, ra) are same as ones of f(ipi(q), ra), Eq. Indeed, the amplitude of the modulations is given by the coefficient 
i?i (see Tab. 0. The correlated character of phonon states is thus related to the feature of the onsite excitation Bloch 
wave ipi(q) and more precisely to the fact that the onsite displacement operator generates an overlap between the 
onsite groundstate 0cm and the first excited state <j>\j,. The agreement between the long range behavior of phonon 
states and ipiil) holds in the harmonic lattice provided that C is not too large. This can be checked by comparing 
the formula /($, ra) = cos(q ra) for a single phonon state in the harmonic chain with the equation Ea ll2l noting that 
for a quadratic form of V, we have D(l, 0) = l/y/2 and thus i?i = 1. 

The agreement between the coefficient R a and the long range modulations of /($,n) is not strictly proved for 
biphonon and triphonon (see Tab[I]and Figs[7Jb-c) and|S£b-c)) because the corresponding branches are too close from 
the unbound phonon bands. Nevertheless we estimate that the agreement roughly holds , whereas in Figs^b) it does 
not because the biphonon gap closes for small q. Provided the a th discrete branch is isolated in the energy spectrum, 
there is no conceptual difference between a phonon bound states and single phonon states because they all come from 
the raising of the translational degeneracy of the onsite eigenstates 4> a ,i- Actually, we conclude that the coefficient R a 
is a good indicator of long range correlations for the a phonon bound states (including phonon states with a = 1) on 
the condition that the hybridization with other phonon states is negligible. Therefore, when R a = 0, the a phonon 
bound states exhibit a finite correlation length provided the a th discrete branch is isolated in the energy spectrum. 
In Fig[7|[b), the space correlation function drops to zero for the large ra and for all q in agreement with R2 equals zero 
(see Tab©. 

As shown in TabQ] R a = 2D(a,0) 2 decreases to zero with increasing a. Above a = 4, R a may be consider as 
physically negligible. By extension of our results about space correlations of the a phonon bound states with a < 4, 
it is expected that the a phonon bound states with a > 4 have a finite correlation length provided their energy 
branch remains isolated from the linear combination bands of the lower energy states, i.e., from the unbound phonon 
states as well as from the combinations of the lower order multi-phonon bound states (for instance, see {21} in Fig[SJ). 
For the specific case ^3 = and A4 > 0, the even values of a verify R a = 0, so the corresponding a phonon bound 
states have also a finite correlation length when their branch is isolated in the spectrum. For a closing biphonon gap, 
a pseudogap was found to persist at the edge of the lattice Brillouin zone fSec lIII|) . In FigHJb), it is shown that the 
biphonon states have a finite length scale in q-range where the pseudogap opens. Indeed, due to hybridization with 
the 2 phonon unbound states, the biphonons at q w have a non-vanishing space correlation function which goes 
beyond the finite lattice size. 
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The phonon bound states with a finite correlation length exhibit a singular feature that is worth characterizing 
more precisely. For that purpose, denote by & a (<l) the wave function of the a phonon bound state with momentum q 
and by E a (q) the dimensionless energy of the a th branch. One introduces the time dependent Wannier state W a (t, n) , 
which is constructed from a combination of the a phonon bound states: 

\W a (t,k) >=^=J2 e ' i(qXk+Eai9)Qt) \^(Q) > ■ (13) 

The subscript k indicates the lattice site where is centered the Wannier transform. In Fig llOl the scalar product 
| < ip a (q)\$> a (q) > | 2 is plotted for a = {1,2,3}. Fixing A3 = and C = 0.05, the parameter A4 is increased from 
zero. The scalar product | < ip a (q)\$ a (q) > | is proved to equal 1 when the gaps that surround the a th branch are 
large, or equivalently when A4 is large. Then, one can reasonably argue that the Bloch wave tp a (q) is a very good 
approximate of the a phonon bound state and the corresponding eigen-energy is approached by a perturbative theory: 



< ^ a (q)W\Mq) >= la + (S- l)( 7 o - CD(Q, 0) 2 ) - C(R a cos(q) + D(a, a)D(Q, 0)). (14) 

When A3 = 0, the correlation function f(ip2,n) is zero for all n (since R2 — 0, see TabUJ, so for the biphonon, we 
obtain E%{q) = 72 + (S — 1)70- Then the energy ^2(9) does not depend on q, ^2(9) = E2, and the Wannier state 
|W2(*> k) > can thus be rewritten as follows: 

\W 2 (t, k) >= e-^^.fclWcM. (15) 

It is obvious that such a state is a localized and time periodic excitation. The classical counterparts of such state 
would be the breather solutions for the classical nonlinear discrete lattice, e.g. A. J. Sievers and S. Takenoi^. These 
classical breather solutions have two important features that are first their spatial localization and second their time 
periodicity with a frequency that is out of both the linear classical phonon branch and its resultant harmonic bands 
(for an exact proof see Ref ll3l) . Our proposition about the classical counterparts could be justified by comparing the 
energies of a localized time periodic Wannier state and the semi-classical quantization of the classical breather orbits, 
in same lattice. In the simple case of zero inter-site coupling C = 0, such a comparison has been performed for the 
onsite double- well potential (see Fig^a) and text in Sec[HJ|. The remarkable agreement obtained at C = 0, allows to 
expect that our proposition holds for weak values of C, at least. Extending our sketch to all a phonon bound states, 
the conditions for a breather-like Wannier state in a KG lattice are twofold: first, the branch of the a phonon bound 
states must be isolated by large surrounding gaps and second, R a = 0. Let us detail the breaking of these conditions. 
When the a th branch is isolated but R a > 0, the i/j a (q) are still some good approximates for the a phonon bound 
states (see in FigllOl for the phonon and triphonon). Then the band width of the a th branch does not vanish because 
of the self-tunneling of tp a) i.e., the intersite coupling hybridizes the Bloch wave tp a (q) with itself. The Wannier 
transform (Ea ll3f) of ip a is a combination of states non-coherent in time, since E a (q) depends on q. It is no longer a 
time-periodic solution. After a certain time At c , the Bloch waves ip a do no longer interfere. This At c varies inversely 
to the band width of the a th branch. On the other hand, when R a = but the a th branch is not isolated, the Bloch 
wave ip a is hybridized with unbound phonon states. Then the Wannier state in Ea ll3l is no longer strongly localized 
as the breather- like Wannier state in Eq^J but it shows some spacial extensions. In addition, E a (q) depends on q 
(as shown in Figl^b), for the biphonon). Hence the Wannier transform is not coherent in time and W a (t, k) is not 
time-periodic. 



V. CONCLUSION 



In the present paper, the pairing of phonon states in the nonlinear KG lattice has been considered with numerics. 
First, the biphonon spectral features have been studied and a pseudogap has been shown to characterize lattices where 
nonlinearity is comparable to phonon band width, whatever are the lattice dimension and the type of nonlinearity. 
Second, we have shown how properties of the KG biphonon depend on the nonlinearity. In the specific case of a strong 
(j) nonlinearity, the biphonon states feature a finite correlation length. The corresponding Wannier transforms are 
breather-like excitations that remain coherent for a duration which increases with the strength of the (j) 4 nonlinearity. 
In contrast, when a 3 nonlinearity predominates, the biphonon exhibits a long range correlation and biphonon 
tunneling is no longer negligible, although it remains weaker than for single phonons. According to our results on 
the triphonon, one may expect that the tunneling of the high energy multi-phonon bound states vanishes for any 
type of nonlinearity, provided that their energy branches avoid the unbound phonon bands. Then some high energy 
breather-like excitations could occur. These states should correspond to the semiclassical quantization of the classical 
breather orbits. 
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Table caption 

Tab. Values of R a = 2 < <f> a ,i\Xi\<j>o i > 2 f° r different a and different nonlinear parameters, ^3 and A4 used 
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Figure captions 

Fig. Energy spectrum of a single He atom embedded into a double-well potential. The rank a of each eigenvalue 
7(a) is given on the X-axis, (a) The semi-classical calculation (square symbols, dashed line) is compared to the 
Hamiltonian diagonalization (circle symbols, solid line) onto the truncated Einstein basis (see text below EqJSJ) with 
a cutoff N = 100. The insert shows the double- well plot versus the displacement (in A), (b) The latter method is 
shown to converge to steady eigenvalues when the basis cutoff N increases (Noo = 2000). 

Fig. |21 Eigen-spectrum of a S = 4 sites chain, with a </> 4 onsite potential. Model parameters are given in Ref. UtI 
Comparison between our numerics (circle symbols) and Ref. (diamond symbols) (see Sec. [nj. The Y-axis unit is 
Ml and its zero is the groundstate energy. The X-axis bears the dimcnsionless first Brillouin zone. 

Fig. 01 Eigen-spectrum of a S = 33 sites chain for parameters: (a) C = 0.05, A 3 = and A 4 = 0.2 and (b) 
C = 0.05, A3 = and A4 = 0.02. The inserts show magnifications of the fundamental branch (left) and the overtone 
region (right). The biphonon branch is marked by ({2}) and the two-phonon band by ({11}). The same tags label 
the inserts. The axis units are same as in Fig. 

Fig. El Same as Fig.^but for different nonlinear parameters: (a) A4 = 0.01 and A3 — 0.13 and (b) A4 = 0.01 and 
A 3 = 0.105. 

Fig. Plot of the energy spectrum of a ID chain, composed of S = 19 atoms for: (a) A4 = 0.2,^3 = 0. and (b) 
A 4 = 0.01, A 3 = 0.13, versus the dimensionless coupling C. Tags are explained in the text. 

Fig. El (color online) On the left hand side, energy spectrum of a 2D square lattice in the region of the biphonon 
for C = 0.025, As — and A4 — 0.1 (top) and for same parameters but A4 = 0.025 (bottom). The lattice size is 
S = 13 X 13. On the right hand side, profiles of the left-hand side spectra along the direction [11]. The inserts show 
the magnifications of the phonon branch (left) and the biphonon energy region (right). 

Fig- Q (color online) Plot of the correlation functions /($,n) for 3 eigenstates of a S = 23 chain. Parameters 
are A4 — 0.2, A3 = 0, C = 0.05, versus the dimensionless distance n. The eigenstates are (a) phonon states , 
(b) the biphonon and (c) the triphonon with for each of them different wave vectors q = (circles), q = An/S 
(triangles), q — 107r/5 (diamonds) , q = I6w/S (squares) and q = 22ir/S (triangles). The inserts in (b) and (c) show 
a magnification of the zero y-axis. 

Fig|S| (color online) Same as Fig. □ but for A 4 = 0.01, A 3 = 0.13 and C = 0.05. 

Fig. El (color online) Same as Fig0but for parameters A4 = 0.02, A 3 = 0. and C = 0.05. The eigenstates are (a) 
the phonon states and (b) the biphonon. 

FigUni Plot of the scalar product | < ip a (q)\$ a (q) > | 2 (see Sec.EJ for C = 0.05 and A 3 — and for q — (solid 
lines) and q = 7r (dashed lines), versus the dimensionless parameter A4. The eigenstates are (a) the phonon states, 
(b) the biphonon states, and (c) the triphonon states. 
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FIG. 4: (2004) L. Proville 
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FIG. 5: (2004) L. Proville: See additional figures 
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FIG. 8: (2004) L. Proville 
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